function [sfunc wfunc] = compWavelet(h, iterations, samples)
  % Compute the scaling and wavelet functions using the 
  % the iterative Cascade algorithm
  % Arguments: iterations (int) - default=6
  %            samples (int)    - default=256
  if nargin <3
    samples = 128; % number of samples over interval  
  end
  if nargin <2
    iterations = 6; % number of iterations of the recursion    
  end
  if nargin <1
    % D4 Daubechies scaling coefficients
    h0 = (1+sqrt(3))/(4*sqrt(2));
    h1 = (3+sqrt(3))/(4*sqrt(2));
    h2 = (3-sqrt(3))/(4*sqrt(2));
    h3 = (1-sqrt(3))/(4*sqrt(2));
    h = [h0 h1 h2 h3];
  end

  g = wavelet_coef(h); % Wavelet filter coefficients
  N = length(h);    % interval is [0,N-1]
  
  % Initial scaling function is a unit pulse for one sample period
  total_samples = (N-1)*samples;
  sfunc = ones(1,total_samples) ./ total_samples; 

  % interpolate scaling and wavelet filters to number of samples
  % by inserting zeros in between the actual coefficients 
  zi = zeros(samples-1,N); 
  him=[h; zi]; hi = reshape(him,1,N*samples); % interpolated h
  gim=[g; zi]; gi = reshape(gim,1,N*samples); % interpolated g

  for i=1:iterations
    % could speed up convolution sum by using FFT in real implementation
    sfunc_new = sqrt(2) * conv(hi,sfunc);
    wfunc_new = sqrt(2) * conv(gi,sfunc);
    
    % decimate scaling function by a factor of 2 on each iteration
    sfunc = sfunc_new(1:2:length(sfunc_new));
    wfunc = wfunc_new(1:2:length(wfunc_new));
  end
  
  xaxis = linspace(0,(N-1),length(sfunc)); % rescale for plot

  figure; 
%  yrange = [-0.5,1.5]; 
xrange = [0,3];
%  do_plot(xaxis,sfunc,xrange,yrange,'scaling',6);  
stairs(xaxis,sfunc);

%  figure; yrange = [-1,1.5];  
%  do_plot(xaxis,wfunc,xrange,yrange,'wavelet',6);
  
function do_plot(x,y,xrange,yrange,funcName,i);
  stairs(x,y);
  xlim(xrange); ylim(yrange);% axis equal;
  msg = sprintf('D4 %s function: iteration %d',funcName, i); title(msg);
  
function g = wavelet_coef(h);
  N=length(h); g = [];
  for i=0:N-1
    gi = (-1)^i * h(N-i);
    g = [g gi];
  end

